#Exposure to Religious Diversity Conditions Apocalyptic Politics – 
#The Failure of the Contact Hypothesis
#Replication File

#Package Install ####
# install.packages("tidyverse")
# install.packages("patchwork")
# install.packages("vtable")
# install.packages("googlesheets4")
# install.packages("interactions")
# install.packages("jtools")
# install.packages("car")
# install.packages("remotes")
# remotes::install_github("ryanburge/socsci")
# install.packages("showtext")
# install.packages("MASS")
# install.packages("ggridges")

#Libraries ####
library(tidyverse)
library(socsci)
library(patchwork)
library(vtable)
library(interactions)
library(jtools)
library(car)
library(MASS)
library(ggridges)
library(showtext)
font_add_google("EB Garamond", "G", regular.wt = 400)
showtext_opts(dpi = 300)
showtext_auto(enable = TRUE)

#Graph theme
theme_results <- function() {
  theme_minimal() %+replace%
    theme(text=element_text(family="G", size=12),
          panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank(),
          legend.position = "bottom",
          plot.title.position = "plot")
}

#Data Import####
resultsdir <- "<insert your own data/results directory here>"
data <- read.csv(paste0(resultsdir, "data.csv"))
                 

#Variable Coding ####
#Religious diversity
data <- data %>% mutate(rdt=q25_1+q25_2+q25_3+q25_4+q25_5)

data <- data %>% mutate(reldiv=(1-(((q25_1/rdt)^2)+((q25_2/rdt)^2)+((q25_3/rdt)^2)+
                                   ((q25_4/rdt)^2)+((q25_5/rdt)^2))))
data <- data %>% mutate(reldivnz=car::recode(reldiv, "0=NA"))

#Herfindahl measure of religious diversity in the state
data <- data %>% mutate(herf2=(((evanrate_2020/100)^2)+((bprtrate_2020/100)^2)+
                               ((mprtrate_2020/100)^2)+((cathrate_2020/100)^2)+((orthrate_2020/100)^2)+
                               ((othrate_2020/100)^2)+((nones_2020/100)^2)))
data <- data %>% mutate(herf2f=1-herf2)

data <- data %>% mutate(pop2020mil=pop2020/1000000)


#Apocalypticism
#Prophecy
data <- data %>% mutate(prophecy=((5-q44_1)+(5-q44_2)+(5-q44_3)+(5-q44_5))/16)

#Evil
data <- data %>% mutate(evil=((5-q42_1)+(5-q42_2)+(5-q42_3))/12)

#End times
data <- data %>% mutate(endtimes=((5-q19_4)+(5-q19_5))/8)

#Christian persecution
data <- data %>% mutate(cp=((4-q16_1)+(4-q16_2)+(4-q16_3)+(4-q16_4))/12)

#Apoc index
data <- data %>% mutate(apoc=(cp+endtimes+evil+prophecy)/4)

#Apoc index minus cp
data <- data %>% mutate(apocno=(endtimes+evil+prophecy)/3)


#Christian nationalism
data <- data %>% mutate(cn=(((5-q14_1)+(5-q14_2)+q14_3+(5-q14_4)+(5-q14_5)+(5-q14_6))/24))


#Demographics
data <- data %>% mutate(female=q74-1)
data <- data %>% mutate(age=2023-q75_1,
                      ed=q80,
                      pid7=q33)
#Partisanship and strength
data <- data %>% mutate(pid7=q33,
                      pidst=car::recode(pid7, "1=4; 2=3;3=4;4=1;5=2;6=3;7=4"))

#Attendance
data <- data %>% mutate(attend5=car::recode(q39, "1:2=5; 3=4; 4=3; 5=2; 6=1"))

#Race
data <- data %>% mutate(white=car::recode(q78_1, "NA=0"),
                      hispanic=car::recode(q78_2, "NA=0"),
                      black=car::recode(q78_3, "NA=0"),
                      asian=car::recode(q78_4, "NA=0"),
                      other=car::recode(q78_5, "NA=0"))

data <- data %>% mutate(race=case_when(white==1 & (black!=1 & hispanic!=1 & asian!=1 & other!=1) ~ 1,
                                     black==1 ~ 2,
                                     asian==1 & (black!=1 & hispanic!=1) ~ 4,
                                     hispanic==1 ~ 3,
                                     other==1 & (black!=1 & hispanic!=1 & asian!=1)~ 5))
data <- data %>% mutate(racef=frcode(race=="1" ~ "White",
                                   race=="2" ~ "Black",
                                   race=="3" ~ "Latino",
                                   race=="4" ~ "Asian",
                                   race=="5" ~ "Other"))


#Need for Chaos
data <- data %>% gather(key="t_remq", value="burn", q26_1, q27_1, na.rm=T)
data <- data %>% gather(key="t_remq", value="remnant", q26_2, q27_2, na.rm=T)

data <- data %>% mutate(chaos2=((5-burn)+(5-remnant))/8)

#Civic skills
data <- data %>% mutate(q30_1=car::recode(q30_1, "NA=0"),
                      q30_2=car::recode(q30_2, "NA=0"),
                      q30_3=car::recode(q30_3, "NA=0"),
                      q30_4=car::recode(q30_4, "NA=0"),
                      skills=q30_1+q30_2+q30_3+q30_4)
                   
#Use of Force
data <- data %>% mutate(force=((5-q21_4)+(5-q21_5))/8)

#Extreme Group Support
data <- data %>% mutate(q32_1c=car::recode(q32_1, "2:3=0"),
                      q32_2c=car::recode(q32_2, "2:3=0"),
                      q32_3c=car::recode(q32_3, "2:3=0"),
                      q32_4c=car::recode(q32_4, "2:3=0"),
                      grprts=(q32_1c+q32_2c+q32_3c+q32_4c))

#Going go Hell
#q38_1, q38_4, q38_3

#Victimhood
data <- data %>% mutate(victim=((5-q18_1)+(5-q18_2)+(5-q18_3)+(5-q18_4)+(5-q18_5)+
                                (5-q18_6)+(5-q18_7)+(5-q18_8))/32)

#RELTRAD ####
data <- data %>% 
  mutate(religpew = q50) %>%  
  mutate(religpew_protestant = car::recode(q51, "14=90; 15=90")) %>% 
  mutate(religpew_baptist = q52) %>% 
  mutate(religpew_methodist = q53) %>% 
  mutate(religpew_nondenom = q54) %>% 
  mutate(religpew_lutheran = q55) %>% 
  mutate(religpew_presby = q56) %>% 
  mutate(religpew_pentecost = q57) %>% 
  mutate(religpew_episcop = q58) %>% 
  mutate(religpew_christian = q59) %>% 
  mutate(religpew_congreg = q60) %>% 
  mutate(religpew_holiness = q61) %>% 
  mutate(religpew_reformed = q62) %>% 
  mutate(religpew_advent = q63)

data <- data %>% 
  mutate(religpew_baptist = car::recode(religpew_baptist,"11=90")) %>% 
  mutate(religpew_methodist = car::recode(religpew_methodist,"6=90")) %>% 
  mutate(religpew_nondenom = car::recode(religpew_nondenom,"6=90")) %>% 
  mutate(religpew_presby = car::recode(religpew_presby,"7=90")) %>% 
  mutate(religpew_pentecost = car::recode(religpew_pentecost,"10=90")) %>% 
  mutate(religpew_episcop = car::recode(religpew_episcop,"5=90")) %>% 
  mutate(religpew_christian = car::recode(religpew_christian,"4=90")) %>% 
  mutate(religpew_congreg = car::recode(religpew_congreg,"4=90")) %>% 
  mutate(religpew_holiness = car::recode(religpew_holiness,"7=90")) %>% 
  mutate(religpew_reformed = car::recode(religpew_reformed,"3=90")) %>% 
  mutate(religpew_advent = car::recode(religpew_advent,"4=90")) 

## Baptist

data <- data %>%
  mutate(sbc = car::recode(religpew_baptist, "1=1; else=0")) %>% 
  mutate(sbc = white + sbc) %>% 
  mutate(sbc = car::recode(sbc, "2=1; else=0"))

data <- data %>%
  mutate(abc = car::recode(religpew_baptist, "2=1; else=0")) %>% 
  mutate(abc = white + abc) %>% 
  mutate(abc = car::recode(abc, "2=1; else=0"))

data <- data %>%
  mutate(ibc = car::recode(religpew_baptist, "5=1; else=0")) 

data <- data %>%
  mutate(bgc = car::recode(religpew_baptist, "6=1; else=0")) 

data <- data %>%
  mutate(mbc = car::recode(religpew_baptist, "7=1; else=0")) %>% 
  mutate(mbc = white + mbc) %>% 
  mutate(mbc = car::recode(mbc, "2=1; else=0"))

data <- data %>%
  mutate(cb = car::recode(religpew_baptist, "8=1; else=0")) 

data <- data %>%
  mutate(fwb = car::recode(religpew_baptist, "9=1; else=0")) 

data <- data %>%
  mutate(gabb = car::recode(religpew_baptist, "10=1; else=0")) 

data <- data %>%
  mutate(obc = car::recode(religpew_baptist, "90=1; else=0")) %>% 
  mutate(obc = white + obc) %>% 
  mutate(obc = car::recode(obc, "2=1; else=0"))

data <- data %>% 
  mutate(evanbap = sbc + abc + ibc + bgc + mbc + cb + fwb + gabb + obc)

##Christian
data <- data %>% 
  mutate(evanchr = car::recode(religpew_christian, "1:3=1; else=0"))

## Methodist
data <- data %>%
  mutate(fmc = car::recode(religpew_methodist, "2=1; else=0")) 

data <- data %>%
  mutate(omc = car::recode(religpew_methodist, "90=1; else=0")) %>% 
  mutate(omc = white + omc) %>% 
  mutate(omc = car::recode(omc, "2=1; else=0"))

data <- data %>% 
  mutate(evanmeth = fmc + omc)

#Reformed
data <- data %>% 
  mutate(evanref=car::recode(religpew_reformed, "2=1; else=0"))

##Non-Denom
data <- data %>% 
  mutate(nd = car::recode(religpew_nondenom, "1:90=1; else=0"))

## Lutheran 
data <- data %>% 
  mutate(mz = car::recode(religpew_lutheran, "2=1; else=0")) %>% 
  mutate(wi = car::recode(religpew_lutheran, "3=1; else=0")) %>% 
  mutate(nalc = car::recode(religpew_lutheran, "5=1; else=0")) %>% 
  mutate(evanluth = mz + wi + nalc)

## Presbyterian

data <- data %>% 
  mutate(evanpres = car::recode(religpew_presby, "2:6=1; else=0"))

## Pentecostal 
data <- data %>% 
  mutate(evanpent = car::recode(religpew_pentecost, "1:90 =1; else=0"))

## Episcopal 
data <- data %>% mutate(acna = car::recode(religpew_episcop, "2=1; else=0"),
                      aoc = car::recode(religpew_episcop, "3=1; else=0"),
                      evanepis = acna + aoc)

## Congregregational
data <- data %>% 
  mutate(evancong = car::recode(religpew_congreg, "2:3=1; else=0"))

## Holiness
data <- data %>% 
  mutate(evanholy = car::recode(religpew_holiness, "1:90 =1; else=0"))

## Advent
data <- data %>% 
  mutate(evanadv = car::recode(religpew_advent, "1:90=1; else=0"))

## None 

## Totaling Up

data <- data %>% 
  mutate(evangelical = evanbap + evanmeth  + evanluth + evanpres + evanpent + 
           evancong + evanholy + evanepis + evanadv + evanchr + evanref) %>% 
  mutate(evangelical = recode(evangelical, "1:4=1; else=0"))

## Making Mainline

data <- data %>% 
  mutate(abc = car::recode(religpew_baptist, "2=1; 4=1; else=0"))

data <- data %>% 
  mutate(epis = car::recode(religpew_episcop, "1=1; else=0"))

data <- data %>% 
  mutate(luth = car::recode(religpew_lutheran, "1=1; 4=1; else=0"))

data <- data %>% 
  mutate(meth = car::recode(religpew_methodist, "1=1; 90=1; else=0"))

data <- data %>% 
  mutate(pres = car::recode(religpew_presby, "1=1; 90=1; else=0"))

data <- data %>% 
  mutate(cong = car::recode(religpew_congreg, "1=1; 3=1; 90=1; else=0"))

data <- data %>% 
  mutate(reform = car::recode(religpew_reformed, "1=1; 3=1; else=0"))

data <- data %>% 
  mutate(mainline = abc + epis + luth + meth + pres + cong  + reform) %>% 
  mutate(mainline = car::recode(mainline, "1:5=1; else=0"))

## Black Protestant 

data <- data %>% 
  mutate(meth = car::recode(religpew_methodist, "3:4=1; else=0"))

data <- data %>%
  mutate(sbc = car::recode(religpew_baptist, "1=1; else=0")) %>% 
  mutate(sbc = black + sbc) %>% 
  mutate(sbc = car::recode(sbc, "2=1; else=0"))

data <- data %>% 
  mutate(nbap = car::recode(religpew_baptist, "3=1; else=0"))

data <- data %>%
  mutate(bpabc = car::recode(religpew_baptist, "2=1; else=0")) %>% 
  mutate(bpabc = black + bpabc) %>% 
  mutate(bpabc = car::recode(bpabc, "2=1; else=0"))

data <- data %>%
  mutate(miss = car::recode(religpew_baptist, "7=1; else=0")) %>% 
  mutate(miss = black + miss) %>% 
  mutate(miss = car::recode(miss, "2=1; else=0"))

data <- data %>%
  mutate(obap = car::recode(religpew_baptist, "90=1; else=0")) %>% 
  mutate(obap = black + obap) %>% 
  mutate(obap = car::recode(obap, "2=1; else=0"))

data <- data %>%
  mutate(ometh = car::recode(religpew_methodist, "90=1; else=0")) %>% 
  mutate(ometh = black + ometh) %>% 
  mutate(ometh = car::recode(ometh, "2=1; else=0"))

data <- data %>% 
  mutate(apos = car::recode(religpew_pentecost, "6=1; 7=1; else=0"))

data <- data %>%
  mutate(open = car::recode(religpew_pentecost, "90=1; else=0")) %>% 
  mutate(open = black + open) %>% 
  mutate(open = car::recode(open, "2=1; else=0"))

data <- data %>%
  mutate(holy = car::recode(religpew_holiness, "90=1; else=0")) %>% 
  mutate(holy = black + holy) %>% 
  mutate(holy = car::recode(holy, "2=1; else=0"))


data <- data %>% 
  mutate(bprot = meth + sbc + nbap + bpabc + miss + obap + ometh + apos + open + holy) %>% 
  mutate(bprot = car::recode(bprot, "1:2=1; else=0"))

## Everything Else

data <- data %>% 
  mutate(catholic = car::recode(religpew, "3=1; else=0"))

data <- data %>% 
  mutate(jewish = car::recode(religpew, "7=1; else=0"))

data <- data %>% 
  mutate(orthodox=car::recode(religpew, "5=1; else=0"))

data <- data %>% 
  mutate(other = car::recode(religpew, "4=1; 8:11=1; 15=1; else=0"))  #maybe also "something else"?

data <- data %>% 
  mutate(none = car::recode(religpew, "12:14=1; else=0"))

data <- data %>% 
  mutate(reltrad = frcode(evangelical == 1 ~ "Evangelical",
                          mainline == 1 ~ "Mainline",
                          bprot == 1 ~ "Black Prot.",
                          catholic == 1 ~ "Catholic",
                          orthodox==1 ~ "Orthodox",
                          jewish == 1 ~ "Jewish",
                          other == 1 ~ "Other Faith",
                          none == 1 ~ "No Religion", 
                          nd == 1 ~ "Non-Denominational",
                          TRUE ~ "Unclassified"))

data <- data %>% mutate(christian=case_when(reltrad=="Evangelical" | 
                                            reltrad=="Mainline" | 
                                            reltrad=="Catholic" | 
                                            reltrad=="Black Prot." | 
                                            reltrad=="Non-Denominational" ~ 1,
                                          TRUE ~ 0))

#Figure 1 – Religious Diversity Acquaintances and their Distribution Across Five Groups ####
data %>% filter(wgt!="NA") %>% mean_ci(rdt, wt=wgt)
data %>% filter(wgt!="NA") %>% mean_ci(reldiv, wt=wgt)

h1 <- data %>% ggplot(aes(x=rdt)) + geom_histogram() +
  theme_results() +
  geom_vline(xintercept=13.5, color="red") +
  labs(x="Total Contacts Listed", y="Count") +
  annotate("text", x=21, y=210, label="Mean=13.5", family="G", size=4) 
h2 <- data %>% ggplot(aes(x=reldiv)) + geom_histogram() +
  theme_results() +
  geom_vline(xintercept=.576, color="red") +
  labs(x="Religious Diversity Measure", y="Count") +
  annotate("text", x=.45, y=210, label="Mean=.58", family="G", size=4) 

h1 + h2
ggsave(file=paste0(resultsdir, "fig1.png"), height=4, width=9)

rm(h1, h2)

#Figure 2 – OLS Estimates of Exposure to Religious Diversity ####
lmrd <- lm(reldiv ~ pop2020mil+herf2f+apoc+attend5+reltrad+age+ed+female+racef+skills, data=data, weight=wgt)
summ(lmrd)
lmrdnz <- lm(reldivnz ~ pop2020mil+herf2f+apoc+attend5+reltrad+age+ed+female+racef+skills, data=data, weight=wgt)
summ(lmrdnz)

plot_summs(lmrd, lmrdnz, inner_ci_level = .90,
           model.names = c("with zeroes", "no zeroes"),
           coefs=c("State population size"="pop2020mil",
                   "State religious diversity"="herf2f",
                   "Apocalypticism"="apoc",
                   "Worship attendance"="attend5",
                   "Mainline Prot"="reltradMainline",
                   "Black Prot" = "reltradBlack Prot.",
                   "Catholic"="reltradCatholic",
                   "Orthdox"="reltradOrthodox",
                   "Jewish"="reltradJewish",
                   "Other faith"="reltradOther Faith",
                   "No religion"="reltradNo Religion",
                   "Non-denominational"="reltradNon-Denominational",
                   "Unclassified"="reltradUnclassified",
                   "Age"="age",
                   "Education"="ed",
                   "Women"="female",
                   "Black"="racefBlack",
                   "Latino"="racefLatino",
                   "Asian"="racefAsian",
                   "Other race"="racefOther",
                   "Civic skills"="skills")) +
  theme_results() +
  geom_vline(xintercept=0, color="red") +
  labs(y="", x="Estimated Effect on Religious Diversity Exposure")
ggsave(file=paste0(resultsdir, "fig2.png"), height=6, width=5)


#Table A2 – OLS Estimates of Religious Diversity Exposure Underlying Figure 2 (p values in parentheses)
export_summs(lmrd, lmrdnz,
             model.names = c("With Zeroes", "Without Zeroes"),
             error_format = "({p.value})", error_pos = "below",
             to.file="docx", file.name = paste0(resultsdir, "FAP1_reldiv.docx"),
             coefs = c("Intercept"="(Intercept)",
                       "State population size"="pop2020mil",
                       "State religious diversity"="herf2f",
                       "Apocalypticism"="apoc",
                       "Worship attendance"="attend5",
                       "Mainline Prot"="reltradMainline",
                       "Black Prot" = "reltradBlack Prot.",
                       "Catholic"="reltradCatholic",
                       "Orthdox"="reltradOrthodox",
                       "Jewish"="reltradJewish",
                       "Other faith"="reltradOther Faith",
                       "No religion"="reltradNo Religion",
                       "Non-denominational"="reltradNon-Denominational",
                       "Unclassified"="reltradUnclassified",
                       "Age"="age",
                       "Education"="ed",
                       "Women"="female",
                       "Black"="racefBlack",
                       "Latino"="racefLatino",
                       "Asian"="racefAsian",
                       "Other race"="racefOther",
                       "Civic skills"="skills"))

rm(lmrd, lmrdnz)

#Figure 3 – Exposure to Religious Diversity Intensifies the Need for Chaos among Apocalyptics ####
lmch <- lm(chaos2 ~ reldiv*apoc+cn+age+ed+female+racef+attend5+reltrad+pidst, data=data, weight=wgt)
summ(lmch)
interact_plot(lmch, pred=apoc, modx=reldiv, interval = T, int.width = .84,
              modx.values = "plus-minus", legend.main = "Religious Diversity") +
  theme_results() +
  labs(x="Apocalypticism", y="Need for Chaos")

ggsave(file=paste0(resultsdir, "fig3.png"), height=4, width=5)


#Table A3 - OLS Estimates for the Need for Chaos Underlying Figure 3 (p values in parentheses)
export_summs(lmch,
             model.names = c("Need for Chaos"),
             error_format = "({p.value})", error_pos = "below",
             to.file="docx", file.name = paste0(resultsdir, "FAP2_chaos.docx"),
             coefs = c("Intercept"="(Intercept)",
                       "Religious diversity"="reldiv",
                       "Apocalypticism"="apoc",
                       "Christian nationalism"="cn",
                       "Age"="age",
                       "Education"="ed",
                       "Women"="female",
                       "Black"="racefBlack",
                       "Latino/a"="racefLatino",
                       "Asian"="racefAsian",
                       "Other"="racefOther",
                       "Mainline"="reltradMainline",
                       "Black Prot."="reltradBlack Prot.",
                       "Catholic"="reltradCatholic",
                       "Orthodox"="reltradOrthodox",
                       "Jewish"="reltradJewish",
                       "Other faith"="reltradOther Faith",
                       "No religion"="reltradNo Religion",
                       "Non-denominational"="reltradNon-Denominational",
                       "Unclassified"="reltradUnclassified",
                       "Partisan strength"="pidst",
                       "Rel diversity * Apoc"="reldiv:apoc"))

rm(lmch)

#Figure 4 – Religious Diversity Intensifies Support for Extreme Group Politics (negative binomial estimates) ####
lmgr <- glm.nb(grprts ~ reldiv*apoc+cn+age+ed+female+racef+attend5+reltrad+pidst, data=data, weight=wgt)
summ(lmgr)
interact_plot(lmgr, pred=apoc, modx=reldiv, interval = T, int.width = .84,
              modx.values = "plus-minus",
              legend.main = "Religious diversity") +
  theme_results() +
  labs(x="Apocalypticism", y="Extreme Group Politics Support")
ggsave(file=paste0(resultsdir, "fig4.png"), height=4, width=5)

#Table A4 – Negative Binomial Estimates of Extreme Group Support Underlying Figure 4 (p values in parentheses)
export_summs(lmgr,
             model.names = c("Extreme Group Support"),
             error_format = "({p.value})", error_pos = "below",
             to.file="docx", file.name = paste0(resultsdir, "FAP3_extreme.docx"),
             coefs = c("Intercept"="(Intercept)",
                       "Religious diversity"="reldiv",
                       "Apocalypticism"="apoc",
                       "Christian nationalism"="cn",
                       "Age"="age",
                       "Education"="ed",
                       "Women"="female",
                       "Black"="racefBlack",
                       "Latino/a"="racefLatino",
                       "Asian"="racefAsian",
                       "Other"="racefOther",
                       "Mainline"="reltradMainline",
                       "Black Prot."="reltradBlack Prot.",
                       "Catholic"="reltradCatholic",
                       "Orthodox"="reltradOrthodox",
                       "Jewish"="reltradJewish",
                       "Other faith"="reltradOther Faith",
                       "No religion"="reltradNo Religion",
                       "Non-denominational"="reltradNon-Denominational",
                       "Unclassified"="reltradUnclassified",
                       "Partisan strength"="pidst",
                       "Rel diversity * Apoc"="reldiv:apoc"))

rm(lmgr)

#Figure 5 – Religious Diversity (Modestly) Intensifies the Need for the Use of Force to Save America Among Apocalyptics ####
lmf <- lm(force ~ reldiv*apoc+cn+age+ed+female+racef+attend5+reltrad+pidst, data=data, weight=wgt)
summ(lmf)
interact_plot(lmf, pred=apoc, modx=reldiv, interval = T, int.width = .84,
              modx.values = "plus-minus",
              legend.main = "Religious diversity") +
  theme_results() +
  labs(x="Apocalypticism", y="Use of Force Required")
ggsave(file=paste0(resultsdir, "fig5.png"), height=4, width=5)


#Table A5 – OLS Estimates of the Use of Force Underlying Figure 5 (p values in parentheses)
export_summs(lmf,
             model.names = c("Use of Force"),
             error_format = "({p.value})", error_pos = "below",
             to.file="docx", file.name = paste0(resultsdir, "FAP4_force.docx"),
             coefs = c("Intercept"="(Intercept)",
                       "Religious diversity"="reldiv",
                       "Apocalypticism"="apoc",
                       "Christian nationalism"="cn",
                       "Age"="age",
                       "Education"="ed",
                       "Women"="female",
                       "Black"="racefBlack",
                       "Latino/a"="racefLatino",
                       "Asian"="racefAsian",
                       "Other"="racefOther",
                       "Mainline"="reltradMainline",
                       "Black Prot."="reltradBlack Prot.",
                       "Catholic"="reltradCatholic",
                       "Orthodox"="reltradOrthodox",
                       "Jewish"="reltradJewish",
                       "Other faith"="reltradOther Faith",
                       "No religion"="reltradNo Religion",
                       "Non-denominational"="reltradNon-Denominational",
                       "Unclassified"="reltradUnclassified",
                       "Partisan strength"="pidst",
                       "Rel diversity * Apoc"="reldiv:apoc"))

rm(lmf)

#Figure 6 – Estimates of Groups Going to Hell Grow with Apocalyptics’ Exposure to Religious Diversity ####
lmhell1 <- lm(q38_1 ~ reldiv*apoc+cn+age+ed+female+racef+attend5+reltrad+pidst, data=data, weight=wgt)
summ(lmhell1)
ipus <- interact_plot(lmhell1, pred=apoc, modx=reldiv, interval = T, int.width = .84,
                      modx.values = "plus-minus",
                      legend.main = "Religious diversity") +
  theme_results() +
  labs(x="Apocalypticism", y="Percent Going to Hell", subtitle="Americans") +
  expand_limits(y=c(0,70))

lmhell4 <- lm(q38_4 ~ reldiv*apoc+cn+age+ed+female+racef+attend5+reltrad+pidst, data=data, weight=wgt)
summ(lmhell4)
ipch <- interact_plot(lmhell4, pred=apoc, modx=reldiv, interval = T, int.width = .84,
                      modx.values = "plus-minus",
                      legend.main = "Religious diversity") +
  theme_results() +
  labs(x="Apocalypticism", y="Percent Going to Hell", subtitle="Christians") +
  expand_limits(y=c(0,70))

lmhell3 <- lm(q38_3 ~ reldiv*apoc+cn+age+ed+female+racef+attend5+reltrad+pidst, data=data, weight=wgt)
summ(lmhell3)
ipgop <- interact_plot(lmhell3, pred=apoc, modx=reldiv, interval = T, int.width = .84,
                       modx.values = "plus-minus",
                       legend.main = "Religious diversity") +
  theme_results() +
  labs(x="Apocalypticism", y="Percent Going to Hell", subtitle="Republicans") +
  expand_limits(y=c(0,70))

ipus + ipch + ipgop + plot_layout(guides="collect", axes="collect") &
  theme_results()
ggsave(file=paste0(resultsdir, "fig6.png"), height=3.5, width=9)


#Table A6 – OLS Estimates of How Many of these Groups are Destined to Hell Underlying Figure 6 (p values in parentheses)
export_summs(lmhell1, lmhell4, lmhell3,
             model.names = c("Americans", "Christians", "Republicans"),
             error_format = "({p.value})", error_pos = "below",
             to.file="docx", file.name = paste0(resultsdir, "FAP5_hell.docx"),
             coefs = c("Intercept"="(Intercept)",
                       "Religious diversity"="reldiv",
                       "Apocalypticism"="apoc",
                       "Christian nationalism"="cn",
                       "Age"="age",
                       "Education"="ed",
                       "Women"="female",
                       "Black"="racefBlack",
                       "Latino/a"="racefLatino",
                       "Asian"="racefAsian",
                       "Other"="racefOther",
                       "Mainline"="reltradMainline",
                       "Black Prot."="reltradBlack Prot.",
                       "Catholic"="reltradCatholic",
                       "Orthodox"="reltradOrthodox",
                       "Jewish"="reltradJewish",
                       "Other faith"="reltradOther Faith",
                       "No religion"="reltradNo Religion",
                       "Non-denominational"="reltradNon-Denominational",
                       "Unclassified"="reltradUnclassified",
                       "Partisan strength"="pidst",
                       "Rel diversity * Apoc"="reldiv:apoc"))

rm(lmhell1, lmhell3, lmhell4, ipus, ipch, ipgop)

#Figure 7 – Religious Diversity (Modestly) Intensifies the Effect of Apocalypticism on Christian Persecution Beliefs and Victimhood ####
lmcp <- lm(cp ~ reldiv*apocno+cn+age+ed+female+racef+attend5+reltrad+pidst, data=data, weight=wgt)
summ(lmcp)
ipcp <- interact_plot(lmcp, pred=apocno, modx=reldiv, interval = T, int.width = .84,
                      modx.values = "plus-minus",
                      legend.main = "Religious diversity") +
  theme_results() +
  labs(x="Apocalypticism", y="Christian Persecution")

lmvic <- lm(victim ~ reldiv*apoc+cn+age+ed+female+racef+attend5+reltrad+pidst, data=data, weight=wgt)
summ(lmvic)
ipvic <- interact_plot(lmvic, pred=apoc, modx=reldiv, interval = T, int.width = .76,
                       modx.values = "plus-minus",
                       legend.main = "Religious diversity") +
  theme_results() +
  labs(x="Apocalypticism", y="Victimhood")

ipcp + ipvic + plot_layout(guides="collect",  axes="collect") & theme_results()
ggsave(file=paste0(resultsdir, "fig7.png"), height=3.5, width=6)


export_summs(lmcp, lmvic,
             model.names = c("Christian\nPersecution", "Victimhood"),
             error_format = "({p.value})", error_pos = "below",
             to.file="docx", file.name = paste0(resultsdir, "FAP6_cp.docx"),
             coefs = c("Intercept"="(Intercept)",
                       "Religious diversity"="reldiv",
                       "Apocalypticism"="apoc",
                       "Christian nationalism"="cn",
                       "Age"="age",
                       "Education"="ed",
                       "Women"="female",
                       "Black"="racefBlack",
                       "Latino/a"="racefLatino",
                       "Asian"="racefAsian",
                       "Other"="racefOther",
                       "Mainline"="reltradMainline",
                       "Black Prot."="reltradBlack Prot.",
                       "Catholic"="reltradCatholic",
                       "Orthodox"="reltradOrthodox",
                       "Jewish"="reltradJewish",
                       "Other faith"="reltradOther Faith",
                       "No religion"="reltradNo Religion",
                       "Non-denominational"="reltradNon-Denominational",
                       "Unclassified"="reltradUnclassified",
                       "Partisan strength"="pidst",
                       "Rel diversity * Apoc"="reldiv:apoc"))

rm(lmcp, lmvic, ipcp, ipvic)

#Table A1 – Unweighted and Weighted Descriptive Statistics ####

datast <- data %>% dplyr::select(reldiv, pop2020mil, herf2f, apoc, cn, age, ed, female, 
                               racef, attend5, reltrad, pidst, pid7, chaos2, skills, 
                               grprts, force, q38_1, q38_4, q38_3, cp, victim, wgt) 

datast %>% dplyr::select(reldiv, pop2020mil, herf2f, apoc, cn, age, ed, female, 
                        reltrad, racef, attend5, pidst, pid7, chaos2, skills, 
                        grprts, force, cp, victim, wgt) %>% 
  drop_na() %>% 
  st(., summ = c('notNA(x)', 'mean(x)',  'sd(x)', 'min(x)', 'max(x)'),
     out="csv", file=paste0(resultsdir, "summtable_noweights.csv"))

datast %>% dplyr::select(reldiv, pop2020mil, herf2f, apoc, cn, age, ed, female, 
                        reltrad, racef, attend5, pidst, pid7, chaos2, skills, 
                        grprts, force, cp, victim, wgt) %>% 
  drop_na() %>% 
  st(., group.weights="wgt", 
     out="csv", file=paste0(resultsdir, "summtable_weighted.csv"))

rm(datast)

#Figure A1 – Distribution of Apocalypticism in the Sample and Among Christians####
data %>% ggplot(aes(x=apoc, y=1)) + 
  geom_density_ridges() +
  geom_density_ridges(data=filter(data, christian==1), fill="dodgerblue2", alpha=.5) +
  theme_results() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  scale_x_continuous(breaks=seq(0,1,.5)) +
  labs(x="Apocalypticism Scale", y="") +
  annotate("text", x=.1, y=2.5, label="Sample", family="G") +
  annotate("text", x=.35, y=4.2, label="Christians", family="G")
ggsave(file=paste0(resultsdir, "figA1.png"), height=4, width=5)


#Figure A2 – Reestimates of Fig 3 Using the Religious Diversity Measure without Zeroes ####
lmch <- lm(chaos2 ~ reldivnz*apoc+age+ed+female+race+attend5+reltrad+pidst, data=data, weight=wgt)
summ(lmch)
interact_plot(lmch, pred=apoc, modx=reldivnz, interval = T, int.width = .84,
              modx.values = "plus-minus",
              legend.main = "Religious diversity") +
  theme_results() +
  labs(x="Apocalypticism", y="Need for Chaos")
ggsave(file=paste0(resultsdir, "figA2.png"), height=4, width=5)

rm(lmch)

#Figure A3 – Reestimates of Fig 4 Using the Religious Diversity Measure without Zeroes ####
lmgr <- glm.nb(grprts ~ reldivnz*apoc+cn+age+ed+female+racef+attend5+reltrad+pidst, data=data, weight=wgt)
summ(lmgr)
interact_plot(lmgr, pred=apoc, modx=reldivnz, interval = T, int.width = .84,
              modx.values = "plus-minus",
              legend.main = "Religious diversity") +
  theme_results() +
  labs(x="Apocalypticism", y="Extreme Group Politics Support")

ggsave(file=paste0(resultsdir, "figA3.png"), height=4, width=5)

rm(lmgr)

#Figure A4 – The Estimated Effects When Substituting Christian Nationalism for Apocalypticism Depend Mightily on Model Specification ####
lmgr3 <- glm.nb(grprts ~ reldiv*cn+age+ed+female+racef+attend5+reltrad+pidst, data=data, weight=wgt)
summ(lmgr3)
ipgr3 <- interact_plot(lmgr3, pred=cn, modx=reldiv, interval = T, int.width = .84,
                       modx.values = "plus-minus",
                       legend.main = "Religious diversity") +
  theme_results() +
  labs(x="Christian Nationalism", y="Extreme Group Politics Support", subtitle="Without Apoc Control")

lmgr4 <- glm.nb(grprts ~ reldiv*cn+apoc+age+ed+female+racef+attend5+reltrad+pidst, data=data, weight=wgt)
summ(lmgr4)
ipgr4 <- interact_plot(lmgr4, pred=cn, modx=reldiv, interval = T, int.width = .84,
                       modx.values = "plus-minus",
                       legend.main = "Religious diversity") +
  theme_results() +
  labs(x="Christian Nationalism", y="Extreme Group Politics Support", subtitle="With Apoc Control")

ipgr3 + ipgr4 + plot_layout(guides="collect") & theme_results()
ggsave(file=paste0(resultsdir, "figA4.png"), height=3.5, width=7)


rm(lmgr3, lmgr4, ipgr3, ipgr4)